home *** CD-ROM | disk | FTP | other *** search
/ Aminet 33 / Aminet 33 - October 1999.iso / Aminet / misc / math / TCalcStats2c.lha / TCalcStats2c / AREXX / Wilcoxon_Sign_Rank.rexx < prev   
Encoding:
OS/2 REXX Batch file  |  1998-08-02  |  8.4 KB  |  423 lines

  1. /* Wilcoxen Sign Rank Test */
  2.  
  3. options results
  4. if ~show('P','TCALC') then do
  5.     address command 'run turbocalc:turbocalc'
  6.     address command 'waitforport TCALC'
  7.     loadflag=1
  8. end
  9. address 'TCALC'
  10. 'DEFPUBSCREEN()'
  11. /* Add-in Rexx Math Library needed for some routines */
  12. signal on syntax
  13. if ~show('l','rexxmathlib.library') then
  14.    call addlib('rexxmathlib.library',0,-30)
  15. if ~show('l','rexxreqtools.library') then
  16.    call addlib('rexxreqtools.library',0,-30)
  17. if ~show('l','rexxsupport.library') then
  18.    call addlib('rexxsupport.library',0,-30)
  19.   /* add to library list */
  20. signal off syntax
  21.  
  22. /* Start Main Routine */
  23. if loadflag=1 then 'Load()'
  24. 'ActivateWindow()'
  25. range=rtgetstring(,"Enter Cell Range for Input","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  26. colon=pos(":",range)
  27. if colon=0 then do
  28.     'Message "Please select a range before executing this script"'
  29.     'DEFPUBSCREEN("Workbench")'
  30.     exit
  31. end
  32.  
  33. /* Find cell references and cell, column numbers */
  34. start_cell=substr(range,1,colon-1)
  35. end_cell=substr(range,colon+1)
  36. start_row=cellrow(start_cell)
  37. end_row=cellrow(end_cell)
  38. start_col=cellcol(start_cell)
  39. end_col=cellcol(end_cell)
  40. NRows=end_row-start_row+1
  41. NCols=end_col-start_col+1
  42. if NCols>2 then do
  43.     'Message "Only Two Columns Allowed!"'
  44.     'DEFPUBSCREEN("Workbench")'
  45.     exit
  46. end
  47.  
  48. /* Get cell reference for output range */
  49. out_cell=rtgetstring(,"Enter Cell Reference for Output","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  50. if out_cell="" then do
  51.     'DEFPUBSCREEN("Workbench")'
  52.     exit
  53. end
  54. if length(out_cell)<2 | datatype(left(out_cell,1),'n') then do
  55.     'Message "Invalid cell reference"'
  56.     'DEFPUBSCREEN("Workbench")'
  57.     exit
  58. end
  59. /* Suppress Screen Redraw to Speed Things Up */
  60. 'Refresh 0'
  61.  
  62. /* Open a small output window on tcalc screen*/
  63. fo=0
  64. CR='0a'x
  65. DisplayMsg="Calculating...Please Wait."||CR||"User input is disabled during calculation."||CR
  66. if open(6Info, 'con:100/0/450/80/Progress/SCREEN TCALC', w) then do
  67.      call writeln(6Info, DisplayMsg)
  68.     fo=1
  69. end
  70. else do
  71.     'Message "TCALC Screen not available for Progress messages"'
  72. end
  73. CALL DELAY(150)
  74.  
  75. /* Get cell references for top cell in each column */
  76. 'SelectCell' start_cell
  77. do col=start_col to end_col
  78.     'GetCursorPos'
  79.     top_cell.col=result
  80.     'Column 1'
  81. end
  82.  
  83. /* Get labels for later use on output */
  84. 'SelectCell' start_cell
  85. 'GetValue'
  86. testlabel=result
  87. testlabel=strip(testlabel)
  88. if datatype(testlabel,'n')=1 then do
  89.     labelflag=0
  90.     do x=1 to NCols
  91.     title.x="Column "||x
  92.     end
  93. end
  94. else do
  95.     labelflag=1
  96.     NRows=NRows-1
  97.     do x=1 to NCols
  98.     'GetValue'
  99.     str=result
  100.     title.x=translate(strip(str),"_"," ")
  101.     'Column 1'
  102.     end
  103. end
  104. if fo then call writech(6Info,"Progress...10 ")
  105.  
  106. /* Get data from cell range */
  107. col=start_col
  108. do x=1 to NCols
  109.     'SelectCell' top_cell.col
  110.     if labelflag=1 then 'CursorDown 1'
  111.     do y=1 to NRows
  112.         'GetValue'
  113.         valtest=result
  114.         if datatype(valtest)='NUM' then do
  115.             'GetValue'
  116.             val=result
  117.             data.x.y=val
  118.         end
  119.         'CursorDown 1'
  120.     end
  121.     col=col+1
  122.     val=0
  123. end
  124. if fo then call writech(6Info,"20 ")
  125. /* Calculate Differences */
  126. d.=0
  127. s.=""
  128. count=0
  129. zero=0
  130. Do y=1 to NRows
  131.     d.y=(data.1.y)-(data.2.y)
  132.     diff.y=d.y
  133.     if (d.y)~=0 then count=count+1
  134.     if (d.y)=0 then do
  135.         zero=zero+1
  136.         s.y=""
  137.         end
  138.     if (d.y)>0 then s.y="pos"
  139.     if (d.y)<0 then do
  140.         s.y="neg"
  141.         d.y=abs(d.y)
  142.         end
  143. end
  144. if fo then call writech(6Info,"30 ")
  145. /*Sort Values */
  146. call Sort()
  147. if fo then call writech(6Info,"40 ")
  148. /*Calculate Ranks */
  149.  
  150. rank.=0
  151. N=count
  152. Do y=1+zero to NRows
  153.     z=y+1
  154.     if (d.y)=(d.z) then do
  155.         cnt=1
  156.         beginrank=y
  157.         avrank=y
  158.         endrank=0
  159.         do until endrank~=0
  160.             y=y+1
  161.             avrank=avrank+y
  162.             z=y+1
  163.             cnt=cnt+1
  164.             if (d.y)~=(d.z) then endrank=y
  165.         end
  166.         avrank=avrank/cnt
  167.         do y=beginrank to endrank
  168.             rank.y=trunc(avrank,2)-zero
  169.         end
  170.     y=y-1
  171.     end
  172.     else rank.y=y-zero
  173.  
  174. end
  175. if fo then call writech(6Info,"60 ")
  176. /* Calculate Rank Sign Sums */
  177. sumpos=0
  178. sumneg=0
  179. Do y=1+zero to NRows
  180.     If s.y="neg" then sumneg=sumneg+(rank.y)
  181.     If s.y="pos" then sumpos=sumpos+(rank.y)
  182. end
  183. T=0
  184. If abs(sumneg)<sumpos then 
  185.     T=abs(sumneg)
  186. else 
  187.     T=sumpos
  188. if fo then call writech(6Info,"90 ")
  189. /* Calculate z */
  190. N=count
  191. aa=T-(N*(N+1)/4)
  192. sdev=sqrt((N*(N+1)*(2*N+1))/24)
  193. zz=aa/sdev
  194. if fo then call writech(6Info,"90 ")
  195. /* Calculate Probability */
  196.  
  197. select
  198.     when zz<(-6.5) then P=0
  199.     when zz>6.5 then P=1
  200.     otherwise
  201.         P=abs(calcp(zz))
  202. end
  203. if fo then do 
  204.     call writeln(6Info,"100")
  205.     call writeln(6Info,"Writing output to window...")
  206. end
  207. /* Output */
  208. dank=0
  209. 'SelectCell' out_cell
  210. 'Put "Wilcoxon Matched-Pairs Signed Ranks Test"'
  211. 'CursorDown 2'
  212. 'GetCursorPos'
  213. top_cell=result
  214. 'Alignment 2'
  215. 'Put "(X-Y)"'
  216. 'CursorDown 1'
  217. Do y=1 to NRows
  218.     'Put' diff.y
  219.     'CursorDown 1'
  220. end
  221. 'CursorDown 1'
  222. 'Put "Sum of Neg. Ranks"'
  223. 'CursorDown 1'
  224. 'Put "Sum of Pos. Ranks"'
  225. 'CursorDown 1'
  226. 'Put "Wilcoxon T:"'
  227. 'CursorDown 1'
  228. 'Put "Ranked Obs. (N):"'
  229. 'CursorDown 1'
  230. 'Put "z:"'
  231. 'CursorDown 1'
  232. 'Put "Prob. (Z<=z):"'
  233. 'SelectCell' top_cell
  234. 'Column 1'
  235. 'ColumnWidth 18'
  236. 'Put "Sorted ABS Values"'
  237. 'CursorDown 1'
  238. Do y=1 to NRows
  239.     'Put' d.y
  240.     'CursorDown 1'
  241. end
  242. 'SelectCell' top_cell
  243. 'Column 2'
  244. 'Alignment 2'
  245. 'Put "Rank"'
  246. 'CursorDown 1'
  247. Do y=1 to NRows
  248.     if (s.y)="neg" then rank.y=-(rank.y)
  249.     if (s.y)="" then do
  250.         'Alignment 2'
  251.         'Put "--"'
  252.         'CursorDown 1'
  253.         end
  254.     else do
  255.         'Put' rank.y
  256.         'CursorDown 1'
  257.     end
  258. end
  259. 'CursorDown 1'
  260. 'Put' sumneg
  261. 'CursorDown 1'
  262. 'Put' sumpos
  263. 'CursorDown 1'
  264. 'Put' T
  265. 'CursorDown 1'
  266. 'Put' count
  267. 'CursorDown 1'
  268. 'Put' format(zz,,4)
  269. 'CursorDown 1'
  270. 'Put' format(P,,6)
  271. 'Refresh 1'
  272. 'Refresh 2'
  273. /*'Message' "Finished"*/
  274. /*indicate the main script is finished*/
  275. DisplayMsg="Cleaning up ...."||CR||"Exiting"
  276. result=writeln(6Info, DisplayMsg)
  277. if result~=0 then do
  278.     /*Wait 3 seconds*/
  279.     CALL DELAY(150)
  280.     /* close window*/
  281.     result=close(6Info)
  282. end
  283. 'DEFPUBSCREEN("Workbench")'
  284. exit
  285.  
  286. /* Procedures */
  287.  
  288. cellrow: procedure
  289. do
  290.     parse arg cell
  291.     do charpos=2 to length(cell)
  292.     if datatype(substr(cell,charpos,1),n) then return substr(cell,charpos)
  293.     end
  294.     return 0
  295. end
  296. Return
  297.  
  298. cellcol: procedure
  299. do
  300.     parse arg cell
  301.     labels="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
  302.     cell=upper(cell)
  303.     len=length(cell)
  304.     val=0
  305. do charpos=1 to len
  306.     if datatype(substr(cell,charpos,1),n) then
  307.     do cell=reverse(substr(cell,1,charpos-1))
  308.     do x=1 to length(cell)
  309.     val=(26**(x-1))*pos(substr(cell,x,1),labels)+val
  310.     end
  311.     return val
  312.     end
  313.     end
  314.     return 0
  315. end
  316. Return
  317. /* It is important to put the exposed array at the end of the next line */
  318. Sort: procedure expose NRows d. s.
  319.  
  320. L=(xtoy(2,int(log(NRows)/log(2))))-1
  321.     Do Until L<1
  322.     L=trunc(int(L/2))
  323.     Do J=1 to L
  324.             Do K=J+L To NRows By L
  325.             I=K
  326.             dumdat=d.I
  327.             dumdats=s.I
  328.             Do while I>L
  329.                 y=I-L
  330.                 If d.y ~> dumdat then Leave
  331.                 d.I=d.y
  332.                 s.I=s.y
  333.                 I=I-L
  334.             End
  335.             d.I=dumdat
  336.             s.I=dumdats
  337.             End
  338.         End
  339.     End
  340. Return
  341.  
  342. syntax:
  343.      if arg(1)='FAIL' then do
  344.         'Message "Library is unavailable."'
  345.         'DEFPUBSCREEN("Workbench")'
  346.         exit
  347.         end
  348.     'DEFPUBSCREEN("Workbench")'
  349.     exit
  350.  
  351. calcp:
  352.     arg g
  353.     factk=1
  354.     pp=0
  355.     itts=0
  356.     term=1
  357.     k=0
  358.     do while abs(term)>(exp(-23))
  359.         term=.3989422804*((-1)**k)*(g**k)/(2*k+1)/(2**k)*(g**(k+1))/factk
  360.         pp=pp+term
  361.         k=k+1
  362.         factk=factk*k
  363.     end
  364.     pp=pp+.5
  365.     if (pp<.0000000001) then pp=0
  366.     return pp
  367.  
  368. Format:  procedure
  369.  
  370.      arg number, before, after
  371.      CallLine = SIGL
  372.      if ~datatype(CallLine, 'N') then CallLine = '??'
  373.  
  374.      /* Make sure we have a number as first (required) argument    */
  375.      if ~datatype(number, 'N') then do
  376.         if number = '' then
  377.            fc = 17     /* Wrong number of arguments           */
  378.         else
  379.            fc = 47     /* Arithmetic conversion error             */
  380.         signal FormatSyntaxError
  381.      end
  382.      num = number + 0
  383.      if before = '' & after = '' then
  384.         return num
  385.      else do
  386.         parse var num integer '.' fraction
  387.         if before = '' then before = length(integer)
  388.         if after = '' then after = length(fraction)
  389.         if ~datatype(before, N) | ~datatype(after, N) then
  390.            do fc = 18
  391.            signal FormatSyntaxError
  392.        end
  393.         if before < length(integer) then do
  394.            fc = 18
  395.            signal FormatSyntaxError
  396.         end
  397.         if after ~= length(fraction) then do
  398.            fraction = trunc(('.'fraction'0') + ('.'copies('0', after)'5'), after)
  399.         if integer<1&integer>-1 then integer=integer
  400.            else integer = integer + (fraction % 1)
  401.            fraction = substr(fraction, 3)
  402.         end
  403.         if fraction >= 0 then
  404.            return right(integer, before)'.'fraction
  405.         else
  406.            return right(integer, before)
  407.      end
  408.  
  409.  FormatSyntaxError:
  410.         if show('F', STDERR) then
  411.            call writeln(STDERR, '+++ Error' fc 'in line' CallLine':' errortext(fc))
  412.         else
  413.            mess='+++ Error' fc 'in line' CallLine':' errortext(fc)
  414.         'Message' mess
  415.         parse source Func .
  416.         if Func = 'FUNCTION' then do
  417.         'DEFPUBSCREEN("Workbench")'
  418.            exit "Err"
  419.         end
  420.         else do
  421.         'DEFPUBSCREEN("Workbench")'
  422.            exit 10
  423.         end